
data {
  int<lower=0> N;         // number of data points
  int<lower=0> N_municipalities;         // number of data points
  int<lower=0> NbDeath[N];
  int<lower=0> Population[N];
  int<lower=0> Municipality_code[N];
  int<lower=0> covid[N]; // standard error of effect estimates
}
parameters {
  real<upper=0> log_hc[N_municipalities];
  real<lower=0, upper=1> hcp[N_municipalities];
  real mu;
  real<lower=0> sigma;
  real<lower=0> lambda;
}
  transformed parameters {
    real hc[N_municipalities];
    hc = exp(log_hc);
  }
model {
  real lambda_eff[N];
  //real hc[N_municipalities];

  mu ~ normal(-5, 5);
  sigma ~ normal(1, 20);
  
  log_hc ~ normal(mu, sigma);
  hcp ~ exponential(lambda);
  lambda ~ gamma(0.001, 0.001);

  for (i in 1:N) {
    lambda_eff[i] = Population[i]*(hc[Municipality_code[i]] + covid[i]*hcp[Municipality_code[i]]);
  }

  //print(mu, mu_p, sigma, sigma_p);
  //print(log_hc);
  //print(log_hcp);
  //print(lambda_eff);

  NbDeath ~ poisson(lambda_eff);
}
generated quantities {
//  real mu_prior;
//  real sigma_prior;
//  real log_hc_prior;
//  real lambda_prior;
//  real hcp_prior;
  int<lower=0> NbDeath_pred[N];
  real lambda_eff[N];


//  mu_prior = normal_rng(-5, 5);
//  sigma_prior = fabs(normal_rng(1, 20));
//  log_hc_prior = normal_rng(mu_prior, sigma_prior);
//  lambda_prior = gamma_rng(0.001, 0.001);
//  hcp_prior = exponential_rng(lambda_prior);
  
  for (i in 1:N) {
    lambda_eff[i] = Population[i]*(hc[Municipality_code[i]] + covid[i]*hcp[Municipality_code[i]]);
  }
  
  NbDeath_pred = poisson_rng(lambda_eff);
}
